#include <iostream>
#include <vector>
#include <string>
#include <fstream>
#include <map>
#include <stdlib.h>
#include <lineages.h>
#include <algorithm>
#include <stdio.h>
#include <sstream>
#include <math.h>
#include <time.h>
#include <shared.h>

using namespace std;

int fitnessType;
int robustType;
int disturbType;

double baseDelRate = 1.0;

int main(int argc, char *argv[])
{
    string realLineagesFile = argv[1];
    fitnessType = atoi(argv[2]);
    robustType = atoi(argv[3]); // must be "3" to use the weights defined in "weightFile"
    disturbType = atoi(argv[4]);
    int simStart = atoi(argv[5]);
    int simLoop = atoi(argv[6]);
    int randSeed = atoi(argv[7]);
    string simFile = argv[8];
    string weightFile = argv[9];
    srand(randSeed);

    double realComplexity,realRobustness;
    double tmpDbl = 0;
    vector<int> termTypes;
    vector<string> termLineages;
    vector<vector<string> > lineagesPool;
    vector<vector<int> > typesPool;

    loadRealLineages(realLineagesFile,termTypes,termLineages);
    Lineages realLin(termTypes,termLineages);

    // load weights
    // user have to make sure the number of weights match the number of types
    ifstream ifs;
    ifs.open(weightFile.c_str());
    int tmpInt = 0;
    while(!ifs.eof())
    {
        tmpInt ++;
        ifs >> tmpDbl;
        ifs.ignore(1000,'\n');
        realLin.types2Weight[0-tmpInt] = tmpDbl;
    }
    tmpInt--;
    if(tmpInt != lin2typeCnt(realLineagesFile))
    {
        cerr << "Expecting " << lin2typeCnt(realLineagesFile) << " Weights" << endl;
        cerr << "Getting" << tmpInt << endl;
        exit(4);
    }

    realLin.getDepthStat();
    realComplexity = realLin.complexity();
    realLin.ruleDeletionRate = baseDelRate / (realLin.termTypes.size() - 1);
    // Here the denominator is (realLin.termTypes.size() - 1)
    // This is correct!! See explanations below, when calculating this metric for the simulated trees.
    realRobustness = realLin.robustness(disturbType,robustType);

    if (disturbType == 2)
    {
        realRobustness = (1 - baseDelRate ) + baseDelRate * realRobustness;
	}
    cout << "Parameters: FitnessType(" << fitnessType;
    cout << ") ; RobustType(" << robustType << ")";
    cout << endl;

    cout << "Lineage Name\tComplexity\tRobustness\tMax depth\tMean depth\tS.D. depth\tfitness" << endl;
    cout << realLineagesFile << "\t" << realComplexity << "\t" << realRobustness << "\t" << realLin.depthMax << "\t";
    cout << realLin.depthMean << "\t" << realLin.depthSD << "\t" << fitness(fitnessType,realRobustness,realComplexity) << endl;

    for(int i=simStart; i<simStart+simLoop; i++)
    {
        lineagesPool.clear();
        typesPool.clear();
        loadLineageCollection(simFile,i,1,typesPool,lineagesPool);

        termTypes = typesPool[0];
        termLineages = lineagesPool[0];

        Lineages simLin(termTypes,termLineages);
        simLin.types2Weight = realLin.types2Weight;
        double simComplexity,simRobustness;
        simLin.getDepthStat();
        simComplexity = simLin.complexity();
        simLin.ruleDeletionRate = baseDelRate / (simLin.termTypes.size() - 1);
        // Here the denominator is (simLin.termTypes.size() - 1)
        // This is correct!! Don't use simLin.internalTypeCnt!!
        // Because rule deletion reflex genetic mutation rate. And in a lineage with more rules, it should have higher chances of mutation/rule deletion.
        // If simLin.internalTypeCnt is used here, a lineage with more rules with have lower mutation rate. That is not what we want.
        // On the other hand, this could of couse be set to some different value, as long as they are reasonable and the same across different simulated lineages based on the same real lineage (Stay the same in the group that with are comparing within) .
        simRobustness = simLin.robustness(disturbType,robustType);

        if(disturbType == 2) // This is obsolete. Don't use it.
        {
            tmpDbl = simLin.internalTypeCnt / realLin.internalTypeCnt ;
            if(baseDelRate * tmpDbl > 1)
            {
                cerr << "Basal Deletion Rate is too high for rule deletion" << endl;
                cerr << "Disturb type (2) is only intended for rule deletion rate < 1" << endl;
                cerr << "Here the effective rule deletion rate for this simulated lineages " << i << " is " << baseDelRate * tmpDbl << endl;
                cerr << "Please consider reduce it" << endl;
            }
            simRobustness = (1 - baseDelRate * tmpDbl) + baseDelRate * tmpDbl * simRobustness;
        }
        cout << i << "\t" << simComplexity << "\t" << simRobustness << "\t" << simLin.depthMax << "\t";
        cout << simLin.depthMean << "\t" << simLin.depthSD << "\t" << fitness(fitnessType,simRobustness,simComplexity)<< endl;
    }

    return 0;
}

